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ABSTRACT 



Context. Recently, new solar model atmospheres have been developed to replace classical ID UTE hydrostatic models and used to for 
example derive the solar chemical composition. 

Aims. We aim to test various models against key observational constraints. In particular, a 3D model used to derive the solar abun- 
dances, a 3D MHD model (with an imposed 10 mT vertical magnetic field), ID NUTE and UTE models from the PHOENIX project, 
the ID MARCS model, and the ID semi-empirical model of Holweger & Miiller. 

Methods. We confront the models with observational diagnostics of the temperature profile: continuum centre-to-limb variations, ab- 
solute continuum fluxes, and the wings of hydrogen lines. We also test the 3D models for the intensity distribution of the granulation 
and spectral line shapes. 

Results. The predictions from the 3D model are in excellent agreement with the continuum centre-to-limb observations, performing 
even better than the Holweger & Miiller model (constructed largely to fulfil such observations). The predictions of the ID theoretical 
models are worse, given their steeper temperature gradients. For the continuum fluxes, predictions for most models agree well with 
the observations. No model fits all hydrogen lines perfectly, but again the 3D model comes ahead. The 3D model also reproduces the 
observed continuum intensity fluctuations and spectral line shapes very well. 

Conclusions. The excellent agreement of the 3D model with the observables reinforces the view that its temperature structure is 
realistic. It outperforms the MHD simulation in all diagnostics, implying that recent claims for revised abundances based on MHD 
modelling are premature. Several weaknesses in the ID hydrostatic models (theoretical and semi-empirical) are exposed. The differ- 
ences between the PHOENIX LTE and NLTE models are small. We conclude that the 3D hydrodynamical model is superior to any of 
the tested ID models, which gives further confidence in the solar abundance analyses based on it. 
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1. Introduction 

Solar model atmospheres are a cornerstone in stellar astronomy. 
The Sun is a natural reference when studying other stars, and 
realistic solar photosphere models are essential to infer solar pa- 
rameters such as its chemical composition. The wealth of solar 
data available can be used to rigorously test and constrain pho- 
tosphere models. This testing provides invaluable insight about 
the model physics and its degree of realism, paving the way for 
building realistic models of other stars. 

With the significant increases in computational power of re- 
cent years, the classical approximations used when modelling 
stellar atmospheres have started to be challenged. Of these ap- 
proximations, the most significant are the assumption of a static 
ID atmosphere with a mixing length type treatment of convec- 
tion, and the assumption of local thermodynamical equilibrium 
(LTE). Although at present no model of a stellar atmosphere 
is able to relax these two assumptions simultaneously, efforts 
have been made to tackle each of these approximations individ- 
ually. On the geometry/convection side, realistic 3D hydrody- 
namical time-dependent simulations of convection have been de- 
veloped and used as models of the solar photosphere (e.g. Stein 
& Nordlund 1998; Asplund et al. 2000a; Freytag et al. 2002; 
Vogler et al. 2004; Carlsson et al. 2004; Caffau et al. 2008b). 
On the radiative transfer side. Short & Hauschildt (2005) have 
computed a ID non-LTE (NLTE) hydrostatic solar model atmo- 



sphere with the PHOENIX code (Hauschildt & Baron 1999) fol- 
lowing the pioneering work in this regard by Anderson (1989). 

The application of 3D solar models to abundance analysis 
(e.g. Asplund et al. 2004; Caffau et al. 2010) has resulted in a 
revised solar photospheric metallicity of Z = 0.0134 (Asplund 
et al. 2009), substantially smaller than previous canonical values 
(e.g. Z = 0.0201 in Anders & Grevesse 1989 and Z = 0.0169 
in Grevesse & Sauval 1998). The realistic treatment of convec- 
tion and velocity fields in the 3D models resulted in an excellent 
agreement between predicted and observed line shapes and bi- 
sectors, not possible with the ID models even with the free pa- 
rameters of micro- and macro-turbulence (Asplund et al. 2000a). 
This agreement is a strong indicator of how realistic the model 
is. Additionally, when compared with observations of the so- 
lar granulation at high spatial resolution, the 3D models cor- 
rectly predict the characteristic size and lifetimes of the granules 
(e.g. Stein & Nordlund 1998; Nordlund et al. 2009). However, 
the results are still controversial because by using a lower solar 
metallicity the previous excellent agreement between solar inte- 
rior models and helioseismology deteriorates significantly (e.g. 
Bahcall et al. 2005; Basu & Antia 2008; Serenelli et al. 2009). 

A criticism of the 3D solar models sometimes raised is that 
while they have been tested against many spectral lines, they 
lack a thorough testing of temperature structure, such as the con- 
tinuum centre-to-limb variation and absolute continuum fluxes 
(Basu & Antia 2008). A correct temperature structure is of the 
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utmost importance to abundance studies. Using a ID horizontal 
and temporal average of the 3D model of Asplund et al. (2000a), 
Ayres et al. (2006) suggest that the 3D model fails to describe the 
observed centre-to-limb variation and its temperature gradient is 
too steep. The first claim is partly dismissed by Koesterke et al. 
(2008), showing that the ID average is not a valid approximation 
of the full 3D model for temperature profiling and that the per- 
formance of the 3D model used by Asplund et al. (2000a) and 
Asplund et al. (2005) in the continuum centre-to-limb variation 
is comparable to that of theoretical ID models - a view corrob- 
orated by Pereira et al. (2008) and Trujillo Bueno & Shchukina 
(2009). However, Koesterke et al. (2008) agree with Ayres et al. 
(2006) in that the temperature gradient of this particular older 
3D model is slightly too steep in the continuum forming layers. 
As mentioned by Asplund et al. (2009), an improved radiative 
transfer treatment resulted in a new 3D model with a more real- 
istic temperature gradient (less steep than before). 

The aim of this work is to systematically test the temper- 
ature structure of several models of the solar photosphere, in- 
cluding 3D and ID NLTE models. We use the new 3D hydrody- 
namical model employed by Asplund et al. (2009), the 10 mT 
3D MHD model of Thaler et al. (in preparation) and recent ID 
NLTE and LTE models from the PHOENIX project (Hauschildt 
& Baron 1999). We also employ the widely used ID MARCS 
model (Gustafsson et al. 2008) and the semi-empirical ID model 
of Holweger & Muller (1974). 

The models used are described in more detail in Sect. 2. To 
compare their temperature structure with the observations we 
use three classical tests: the continuum centre-to-limb variation 
in Sect. 3, the absolute continuum flux distribution in Sect. 4 and 
the wings of hydrogen lines in Sect. 5. In addition, we also test 
the 3D model against observations of the continuum intensity 
distribution and A/n^j in Sect. 6. Sect. 7 contains a comparison of 
the predicted and observed line shapes, including inferred abun- 
dances, bisectors and line shifts, for a sample of Fe i and Fe ii 
lines. Our conclusions are given in Sect. 8. 



2. Model atmospheres 

2.1. 3D Model 

We use the same 3D hydrodynamical model solar atmosphere 
that was adopted by Asplund et al. (2009) for the derivation 
of photospheric solar abundances. The solar surface convection 
simulation was performed using the 3D, radiative, hydrodynam- 
ical, conservative, stagger-code (Nordlund & Galsgaard 1995). 
In the simulation, the equations for the conservation of mass, 
momentum, and energy are solved together with the radiative 
transfer equation for a representative volume of solar surface 
(6x6x3.8 Mm^) on a Cartesian mesh with 240^ numerical reso- 
lution. The horizontal grid is equidistant, while the vertical depth 
scale has a non-constant spacing optimised to better resolve the 
layers at the photospheric transition where temperature gradients 
are at their steepest. Open, transmitting, boundaries are assumed 
at the top and bottom of the simulation domain, and periodic 
boundary conditions are enforced horizontally. It is important for 
the lower boundary to be transmitting, to avoid homogenising 
the otherwise highly asymmetric (between up and down) con- 
vective flows. 

The simulation domain completely covers the Rosseland op- 
tical depth range -5 < logTRoss ^ 7. The radiative transfer equa- 
tion is solved using a long-characteristics Feautrier-like scheme 
down to rRoss~300, and the diffusion approximation is employed 
in the deeper layers. The main improvement over the simulation 



of Asplund et al. (2000a) is the treatment of opacities and line- 
blocking in particular. In both the old (2005) and new (2009) 
3D simulations, continuous and line opacities are included via a 
statistical method, called opacity binning or multi-group method 
(Nordlund 1982): wavelengths are sorted into opacity bins ac- 
cording to the strength of the opacity and the corresponding LTE 
source functions are added together within each bin. In the orig- 
inal binning scheme, Nordlund (1982) used the Rosseland aver- 
age K() of the opacities in the continuum bin, scaled by a con- 
stant factor for each of the other bins, /q - /cqIO-'*^, in prac- 
tice assuming Ax = 1 and /' = 0, . . ., 11. The transition to free 
streaming for optical depths in the continuum bin to<s1 is en- 
sured by an exponential bridging (in tq ) to an intensity-weighted 
mean opacity. The multi -group method was further developed by 
Skartlien (2000), who relaxed the approximation of kj just be- 
ing a scaled kq, and instead computed the actual Rosseland av- 
erage for each bin. This also ensures that the Rosseland mean of 
the bin-wise opacities converges to the actual Rosseland mean 
of the monochromatic opacities. The new 3D model has fur- 
ther been improved by sorting opacities into bins not only ac- 
cording to opacity strength but also according to wavelength, 
and allowing arbitrary bin sizes; a similar binning criterion has 
been implemented by Caff'au et al. (2008a) in the CO^BOLD 
code. For the present simulation, continuous opacities are taken 
from Gustafsson et al. (1975, and subsequent updates) and line 
opacities from the latest MARCS stellar atmosphere package 
(Gustafsson et al. 2008). The solar chemical composition by 
Asplund et al. (2005) and the equation-of-state by Mihalas et al. 
(1988, and subsequent updates) are adopted. 

The positions of the bin borders are then optimised with re- 
spect to a monochromatic radiative transfer calculation for the 
simulation's average temperature and density stratification taken 
on surfaces of constant Rosseland optical depth. The generalisa- 
tion of the bins and the optimisation reduce the differences be- 
tween the radiative heating of the monochromatic and the binned 
solution by a factor of five, to within less than 1%. The effect 
of radiative heating and cooling on the simulation is therefore 
faithfully reproduced by this opacity-binning, for the average 
atmosphere. We have also performed a 3D monochromatic ra- 
diative transfer calculation on the full 3D simulation cube for a 
single snapshot, and found a < 5 K (0.08%) difference in T^s 
between the monochromatic and the binned solution. Before be- 
ing subjected to scientific tests, the simulation has been fully 
relaxed, and has subsequently run for the 45 solar minutes used 
here (90 snapshots). The relaxation process included extracting 
energy from radial /7-modes, and ensuring that the total flux is 
statistically constant with depth, that no drifts are present in the 
thermodynamical quantities at the bottom boundary, and that the 
vertical grid's resolution is enough for the radiative transfer in 
the photosphere. 

For the calculations presented here the original simulation 
was interpolated to a coarser 50x50x82 resolution (with a finer 
vertical depth scale) to save computing time. The effective tem- 
perature of the 90 snapshots used is Teff = 5778 + 2 K, close to 
the observed value of 5777 ± 3 K (Willson & Hudson 1988). 

2.2. 3D MHD Model 

In addition to the 3D simulation, we also test a 3D magneto- 
hydrodynamical (MHD) simulation from Thaler et al. (in prepa- 
ration). This simulation was performed using the same code 
and physical ingredients of the 3D hydrodynamical model of 
Asplund et al. (2009). It uses 240x240 grid points horizontally 
and 220 vertically. It has the same horizontal size of 6x6 Mm^ as 
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Fig. 1. Temperature structure of the 3D and ID models, plotted against the optical depth at 500 nm. For the 3D models structure 
represents the temporal and spatial mean (over T500 iso-surfaces). Bottom panels: differences between the 3D model and a given 
model (legend according to the top panels). 



the 3D hydrodynamical model, but a slightly shorter vertical ex- 
tension of 3.34 Mm. It extends 2.7 Mm into the convection zone 
and reaches 0.645 Mm up the photosphere. As in the 3D hydro- 
dynamical model, the horizontal grid is equidistant and the verti- 
cal depth scale optimised to resolve the photosphere. The bound- 
ary conditions are the same as in the Asplund et al. (2009) sim- 
ulation. On a thermodynamically relaxed hydrodynamical snap- 
shot, a vertical magnetic field of 10 mT was overimposed. The 
magnetic field is kept vertical at the bottom boundary and tends 
toward a potential field at the top. This configuration was run 
for 120 min of solar time. For the results presented here, a se- 
quence of 38 min was used, taken after this simulation had run 
for 68 min. 

The choice of 10 mT for B was made as it is a reasonable 
value for the quiet sun's mean field strength (e.g. Trujillo Bueno 
et al. 2006) and because it allows a comparison with the middle 
MHD model of Fabbian et al. (2010, 2012). 



2.3. 1D Models 

We use two types of ID models: theoretical and semi-empirical. 
The semi-empirical ID model of Holweger & Muller (1974) was 
built from a range of observables to reproduce the mean physical 
quantities of the solar photosphere. Most importantly, it was con- 
structed to follow the observed continuum centre-to-limb varia- 
tion between 0.5-300 /im and the line depths of a;900 spectral 
lines. This is an important detail to note when considering our 
centre-to-limb variation comparison, although the observations 
we employ are more recent than the ones available when the 
Holweger & Muller model was built. Historically, the Holweger 
& Muller model has been the atmosphere of choice when de- 
riving solar abundances. It assumes hydrostatic equilibrium and 
does not explicitly include convection. 

Of the ID hydrostatic theoretical model atmospheres we 
include the LTE, line-blanketed solar MARCS model, which 
is a reference for the MARCS grid of model atmospheres 
(Gustafsson et al. 2008). We also include an LTE and an NLTE 
model from the PHOENIX project (Hauschildt et al. 1999). 



These two models have been computed for the solar abundances 
of Asplund et al. (2005) and the same input physics as for 
the PHOENIX Gaia grid (Brott & Hauschildt 2005). They dif- 
fer only in their treatment of atomic level populations with the 
NLTE model having been computed with a NLTE treatment of 
H, He, C, N, O, Mg and Fe. 

To ensure consistency when using our line formation code, 
opacities and equation of state, for the ID models we took the 
T(t) relation from their respective references and integrated Pg^s 
in optical depth assuming hydrostatic equilibrium to obtain the 
pressures and densities that yield the same T(t) when using our 
opacities and equation of state. We note that this is an often over- 
looked source of error when comparing results for supposedly 
the same model atmosphere since this pressure-integration is not 
always carried out. In fact the Holweger & Muller (1974) model 
is essentially only an updated version of the Holweger (1967) 
model with a new pressure-integration due to their effects on 
opacities and thermodynamics. 



2.4. Mean stratification 

In Fig. 1 we compare the mean temperature structure of the 
3D hydrodynamical model with other models. On the left, with 
the 3D MHD model of Thaler et al. (in preparation), the ID 
Holweger & Muller and MARCS models (left) and on the right 
with the old 3D model of Asplund et al. (2000a), and the 
PHOENIX LTE and NLTE models. The mean structure of the 
3D models was calculated by averaging over tsoq iso-surfaces. 

The 3D model of Asplund et al. (2000a) is provided here for 
a quick comparison. This older model was run with a precursor 
of the STAGGER-CODE, and was used to derive the solar abundances 
of Asplund et al. (2005). The largest difference between the old 
and new models is the treatment of radiation. The new model has 
an improved multi-group opacity scheme with 12 bins, whereas 
the old model had a more approximate scheme and used only 4 
bins. From Fig. 1 one can see that the major consequence of the 
improved scheme is a warmer upper photosphere, but unchanged 
at the top of the domain, resulting in a shallower temperature 
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Fig. 2. Continuum centre-to-limb observations in the visible and 
near-infrared, from Pierce & Slaughter (1977), Pierce et al. 
(1977) and Neckel & Labs (1994). In the wavelength region 
where the sets overlap we use only Neckel & Labs (1994) for 
our comparisons. 

gradient. This difference in temperature gradient will be most 
noticeable in the continuum centre-to-limb variations, which we 
discuss below. 

The effect of NLTE in the PHOENIX models seems to be 
a cooling of the outer layers (~ 50 K), with minor differences 
at other depths. This NLTE cooling goes in the opposite direc- 
tion of the NLTE effects of other PHOENIX solar NLTE models 
(Short & Hauschildt 2005, 2009), where the NLTE effects cause 
a warming in the outer layers as naively expected from reduced 
line blanketing and surface cooling. This discrepancy seems to 
be associated with a different choice of atomic species treated in 
NLTE. 

When debating on the advantages and disadvantages of em- 
ploying a full 3D analysis to stellar photospheres, an important 
question to ask is: is the advantage brought by the 3D treatment 
of convection merely a realistic temperature stratification, or are 
the spatial and temporal inhomogeneities essential to derive ac- 
curate abundances? To answer these questions we included an- 
other ID model in the testing, corresponding to the spatial and 
temporal mean structure of the 3D model (shown in Fig. 1). The 
horizontal averaging was done on surfaces of equal optical depth 
rather than geometrical height. This model, hereafter the (3D) 
model, will enable us to disentangle the effects of the mean tem- 
perature structure from the full 3D treatment. 

3. Continuum centre-to-limb variations 

3.1. Context 

The centre-to-limb variations (CLV) of continuum intensities 
provide a sensitive probe of the solar photosphere. Because the 
continuum intensity is proportional to the local source function 
of continuum forming regions, its CLV is a measure of the tem- 
perature variation with depth (the closer to the solar limb, the 
higher up in the atmosphere). The variation of depth can be ex- 
pressed in terms of ju = cos 6, where is the heliocentric viewing 
angle. Normalising /(ju) by the disk-centre value I(fi = 1), one 
has a measure of the temperature gradient of the photosphere 
around the continuum forming layers. This provides a robust test 
of models. 



3.2. Observations 

We make use of the CLV observations of Neckel & Labs (1994) 
and Pierce et al. (1977). They cover respectively the wavelengths 
between 303-1099 nm and 740.4-2401.8 nm, as shown in Fig. 2 
for 0.1 < fi < 0.9. We compare the models with Neckel & 
Labs (1994) for A < 1099 nm and Pierce et al. (1977) for 
A > 1099 nm. Although not used in our comparison, the ob- 
servations of Pierce & Slaughter (1977), covering the range 
303.3-729.7 nm, are also plotted in Fig. 2 and agree very well 
with Neckel & Labs (1994). 

Other CLV observations of longer wavelengths exist, such as 
Spickler et al. (1996). However, as is visible in Fig. 3 of Spickler 
et al. (1996), there is a considerable scatter between different ob- 
servations, especially for A> 4 fim. Because of these uncertain- 
ties we do not include them in this comparison. 

3.3. Results and discussion 

The model predictions were computed using our 3D LTE line 
formation code, which was used to obtain the predicted inten- 
sity at different inclinations. The intensities were computed for 
nine different values of f^ = cos 6 and for each inclination ex- 
cept the vertical they were averaged over four ^-angles in the 
3D case. These intensities were interpolated in jj for the incli- 
nations shown. For the 3D model the intensities were computed 
for each of the 90 snapshots, the final value being the spatial and 
temporal average. 

Results for the visible and near-infrared continuum centre- 
to-limb variation are shown in Fig. 3. A more compact compar- 
ison with Neckel & Labs (1994) is made in Fig. 4, where we 
plot the (normalised) absolute value of the difference between 
models and observations, averaged over wavelength for each yu 
value. Because of uncertainty in the observations (from the dif- 
ficulty in finding continuum regions), this average is limited to 
the 400 - 1099 nm wavelength region. 

The results show an excellent agreement between the 3D 
hydrodynamical model of Asplund et al. (2009) and the ob- 
servations, particularly when comparing with Neckel & Labs 
(1994). This agreement is visibly better even than that of the 
ID Holweger & Miiller model, which is quite remarkable given 
that it was empirically constructed to fit the centre-to-limb vari- 
ations. It should be noted, however, that the observational and 
atomic data have improved much since the construction of the 
Holweger & Miiller model. The application of modern opaci- 
ties to the Holweger & Miiller T(t) stratification therefore re- 
sults in discrepancies not present at its development in 1974. 
We note also that the spectral resolving power of the solar at- 
las employed by Holweger & Miiller was poorer than the at- 
lases available today, leading to less steep line cores and con- 
sequently a too high temperature was inferred from the spectral 
inversion process; rectifying this shortcoming results in a tem- 
perature structure very similar to the mean of the here employed 
3D model (N. Grevesse, 2009, private communication). With the 
infrared observations the agreement with the models is slightly 
worse. However, observations at these wavelengths seem more 
uncertain, as can be seen by the increased scatter between data 
points. It is likely that this additional uncertainty affects the 
agreement with the models. Indicative of this is the region be- 
tween 1500-1850 nm, where the observations show some scatter 
and are consistently higher than the model predictions. 

The inclusion of a 10 mT magnetic field in the 3D model 
somewhat degrades the agreement with continuum centre-to- 
limb variation as seen in Fig. 3 and 4. This difference can be 
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Fig. 3. Centre-to-limb variations in the continuum intensity. Top panels: comparison with the visible/infrared observations of Neckel 
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traced to the shallower temperature gradient (Fig. 1) in the MHD 
model. The 3D MHD model performs slightly worse than the 
Holweger & Muller model. 

The agreement with the theoretical ID models is not as 
good. It is interesting to note in Fig. 4 that LTE models of 



MARCS and PHOENIX have the same trend with ju, although 
the PHOENIX model performs slightly better. The results for 
the PHOENIX NLTE model depart only slightly from the LTE 
model results. The NLTE cooling of the outer layers seen in 
Fig. 1 causes a slightly steeper temperature gradient, which leads 
to a worse agreement with the observed centre-to-limb varia- 
tions. The overall structure and dependence with ji remains es- 
sentially the same for both PHOENIX models as well as for the 
MARCS model, as seen in Figs. 3 and 4, due to the similarity in 
T{t) for -2 < log T < 0, the layers largely tested with continuum 
CLV. 

Compared to other models, the differences between the 3D 
and (3D) models are small, meaning that the mean temper- 
ature gradient is the main driver of the continuum CLV be- 
haviour Nevertheless, the 3D model predictions agree even 
closer with the observations, confirming the results of Koesterke 
et al. (2008), although we find a smaller "3D-(3D)" difference. 
Looking at Fig. 3, the (3D) model lies slightly below the 3D 
model, in other words the effect of the atmospheric inhomo- 
geneities increases I{fi)/I(fi =1). One would therefore expect 
that if spatial and temporal inhomogeneities were added to the 
Holweger & Muller model, its predictions would lie further 
above the observations. This indicates that the temperature gra- 
dient of the Holweger & Muller model is too shallow compared 
to the Sun. 

We compare also with the old 3D model. While this is a re- 
alistic model that reproduces the observed line shifts and shapes 
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(Asplund et al. 2000a), its steeper temperature gradient has a no- 
ticeable effect on the continuum CLV. Its predictions fare worse 
when compared to the observations (but still better than the ID 
models). We do not use this old model in the other observational 
tests. 

In summary, the confrontation with continuum centre-to- 
limb observations reveal that the 3D hydrodynamical model 
is very realistic. The (3D) hydrodynamical model, the 3D 
MHD and the Holweger & Miiller model all perform slightly 
worse with too little limb-darkening, while the MARCS and 
PHOENIX ID theoretical models all predict much too strong 
centre-to-limb variation due to a too steep temperature gradient. 



4. Absolute flux distribution 

4.1. Context 

An independent test of the temperature structure of solar mod- 
els is provided by the observed absolute continuum fluxes. They 
provide an absolute temperature scale for the photospheric lay- 
ers. 
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Fig. 6. Continuum flux differences between the models and the 
observations. For A < 450 nm the observations are not very reli- 
able because of difliculties in continuum placement. The feature 
at A X 950 nm is likely to be caused by uncorrected telluric ab- 
sorption in the observations. 



4.2. Observations 

Several observations of the absolute solar flux/irradiance are 
available, either space-based (e.g. Colina et al. 1996; Thuillier 
et al. 2004) or Earth-based (e.g. Kurucz et al. 1984; Brault & 
Neckel 1987; Kurucz 2005, all obtained with the high-resolution 
NSO/Kitt Peak Fourier transform spectrometer). The space- 
based observations provide a spectrum clean of terrestrial ab- 
sorption features, but at the cost of a lower spectral resolution. 

To extract the observed continuum fluxes we used an abso- 
lute flux atlas divided by a normalised flux atlas (for the points 
deemed to be near the continuum). We used the Kurucz (2005) 
irradiances and normalised flux atlases, available for the interval 
of 300-1000 nm. They are a recent reduction of the data used to 
produce the atlas of Kurucz et al. (1984), and have been carefully 
adjusted to remove the telluric absorption features. The choice of 
using the Kurucz (2005) atlases instead of space-based observa- 
tions was made because of its consistency between absolute and 
normalised fluxes. 

If another irradiance atlas was used, the slight differences in 
line strengths or wavelength mismatches between the irradiance 
and normalised flux atlases would cause some scatter in the con- 
tinuum fluxes, which would have to be smoothed out (see Ayres 
et al. 2006, where such an approach is followed). To avoid these 
uncertainties we use the Kurucz (2005) atlases that, being pro- 
duced from the same data, have a consistent continuum deter- 
mination and wavelength calibration between the irradiance and 
normalised flux atlases. 

The observed continuum flux was obtained as follows. First 
the irradiance is converted to flux at the solar surface using the 
multiplicative factor of [(1 AU)/RQf- - 46202. Then the wave- 
lengths of the (near) continuum points in the normalised flux 
are identified. These are defined as F^ > 0.99 F™, where F™ 
is a local maximum in 5 nm windows. Finally the absolute flux 
is divided by the normalised flux for the continuum high wave- 
lengths. This ratio is linearly interpolated to a coarser resolution 
of 1 nm. A few spurious points were manually removed. The 
Kurucz (2005) irradiance has been normalised to the total solar 
irradiance of Thuillier et al. (2004). Following the discussion in 
Ayres et al. (2006) we have readjusted the continuum fluxes by 
-0.4%, to account for the more accurate total solar irradiance 



of Kopp et al. (2005). In Fig. 5 we show the original observed 
fluxes along with our derived continuum fluxes. 

4.3. Results and discussion 

The predicted continuum fluxes were computed with our LTE 
line formation code. The disk-integrated fluxes are computed us- 
ing eight /i-angles and four (/^-angles, a total of 32 angles; in ad- 
dition the disk-centre intensity was calculated. The //-angles and 
weights are taken from a Gaussian quadrature and are thus not 
identical to those used for the centre-to-limb study presented in 
Sect. 3. For the 3D models the fluxes are computed and spatially 
and temporally averaged for the 90 simulation snapshots. For ID 
models we used the same nine ju-angles, including disk centre. 

The resulting fluxes for all models are shown in Fig. 5. A dif- 
ferential comparison is also shown in Fig. 6. There is an excess 
of flux for A < 370 nm in the observed continuum fluxes when 
compared with the original flux (with lines) or any of the mod- 
els. This difference, also evident in Fig. 6, seems to be caused 
by a too high continuum placement in the normalised Kurucz 
(2005) flux atlas (e.g., higher than in Kurucz et al. 1984; Brault 
& Neckel 1987). Being a region very crowded with lines, the 
discrepancy between observations and models likely arises from 
the difficulty in finding the continuum level (which is systemati- 
cally overestimated), and not from a failure of the models. 

Overall, the ID MARCS model is the best at reproducing the 
observed flux, although the differences between different mod- 
els are small. The Holweger & Muller predictions are consistent 
with a correct Teff, but its flux distribution has a different shape. 
Between 350-450 nm it has less flux, but beyond A x 500 nm 
it shows an excess flux when compared with the observations. 
The 3D model consistently predicts slightly less flux than the 
observations, but nevertheless has a good agreement and repro- 
duces the flux distribution well. The PHOENIX LTE model be- 
haves similarly to the 3D model, but with less flux at the peak 
of the distribution. The differences between the NLTE and LTE 
PHOENIX models are very small, not surprising given that very 
little flux comes from the regions where the NLTE cooling is 
more efficient. 

In the comparison between (3D) and 3D, one finds consid- 
erably less flux coming from the (3D). This is because the T500 
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Fig. 5. Absolute fluxes for the models and observations. Top left: Comparison of the original fluxes of Kurucz (2005), our derived 
continuum fluxes (see text), and the predicted continuum fluxes from the 3D model. Other panels: observed continuum fluxes and 
predicted continuum fluxes for several models. 



iso-surface averaging per construction does not preserve the ef- 
fective temperature. In the 3D model a reasonable amount of flux 
will be emitted from the snapshots with a higher Teff , asF cc T*^. 
The (3D) model is not averaged in T* and will have a lower r^ff . 
This is perhaps one of the strongest arguments against construct- 
ing modified ID models that recover the mean structure of real- 
istic convection: T^ff will not be preserved. One can argue that 
the (3D) model could be averaged in T'^. The issue of how to av- 
erage 3D models will be addressed in a forthcoming paper; tests 
performed so far indicate that for the purposes of line formation 
there is little difference between T- and T^-averaged models. 



Table 1. Root mean square differences between the models and 
the observations, between 375-975 nm. 



Model 


AF'/AF-;^ 


3D Model 


1.00 


3D MHD Model 


0.78 


ID Holweger & Miiller 


3.31 


ID MARCS 


0.67 


ID PHOENIX LTE 


2.96 


ID PHOENIX NLTE 


3.98 


(3D> Model 


15.00 



The flux differences in Fig. 6 are quantified in Table 1, where 
we show the root mean square differences between the mod- 
els and observations, summed over the region of 375-975 nm, 
and normalised to the results from the 3D model. Again one 
can see that the MARCS model gives the best flux predictions, 
followed by the 3D MHD and 3D HD model. The PHOENIX 
and Holweger & Miiller models perform similarly while inter- 
estingly the (3D) model shows the largest differences. 



5. Hydrogen lines 

5.1. Context 

The first two lines in the hydrogen Balmer series lines. Ha 
and HjS, have traditionally been used for temperature determi- 
nation in late-type stellar atmospheres. Their wings are typically 
formed in the region around -2 < log rgoo < 0.5, deeper than 
most of the other spectral lines in late-type stars (Fuhrmann et al. 
1993), while their cores are formed in the chromosphere. The 
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fact that for late-type stars these lines are more sensitive to tem- 
perature than they are to gravity, metallicity or indeed the hydro- 
gen abundance has established them as a popular tool for tem- 
perature determinations (Fuhrmann et al. 1993, 1994; Barklem 
et al. 2002; Behara et al. 2009). The shapes of hydrogen lines, 
because of their large depths of formation, are however depen- 
dent on the convection efficiency (e.g. Ludwig et al. 2009). 

Also relevant to this work are the lines of the Paschen series. 
With /iiow = 12.088 eV, these lines are formed even deeper than 
the Balmer series. However, owing to their longer wavelengths 
and lower level populations (and consequently weaker wings), 
they have not been used as extensively as Ha and HyS. We include 
in our comparisons the Pay and Pa/3 lines, which are just inside 
the wavelength range of the high-resolution FTS atlases of the 
Sun. 

Departures from LTE can be important for hydrogen lines 
(Przybilla & Butler 2004; Barklem 2007). While the far wings 
are formed in deep, hot regions where LTE conditions largely 
prevail, the LTE assumption is no longer valid for the cores of the 
strong lines. However, Barklem (2007, hereafter PB07) shows 
that in the case of Balmer lines, the NLTE effects can even ex- 
tend into the line wings, depending on the adopted rates for col- 
lisions with H atoms. In their LTE and NLTE study of solar H 
lines, Przybilla & Butler (2004) find a weakening of the NLTE 
effects for the Paschen lines. To obtain realistic line profiles, we 
perform NLTE line synthesis as detailed in Sect. 5.3. 

5.2. Observations 

For the solar observations of hydrogen lines we use the high- 
resolution FTS normalised flux atlas of Kurucz (2005), and the 
FTS disk-centre atlas of Brault & Neckel (1987). For our com- 
parison we study the Balmer lines He and H/3 along with the 
Paschen lines Pay and Pa/3. 

In addition to the disk-centre atlas, Brault & Neckel (1987) 
have also made available a flux atlas. For the lines considered 
here it is essentially identical to the Kurucz (2005) atlas, which 
is not surprising since it was produced from the same raw FTS 
data. However, we would like to point out that the flux atlas 
of Brault & Neckel (1987) has a slightly different normalisa- 
tion, with a lower continuum level. This difference amounts to 
» 0.4-0.6% of the normalised flux. For a consistent comparison 
of the flux and disk-centre intensity profiles, we have increased 
the continuum level of the Brault & Neckel (1987) disk-centre 
atlas by 0.5%, so that it matches that of the Kurucz flux atlas. 
This choice of continuum will not affect the differential results 
between models, only their relative standing to the observations. 
Our choice of continuum falls on the Kurucz atlas because it is a 
more recent (and hopefully more accurate) reduction of the same 
data. However, the 0.5% continuum error, which corresponds to 
a Teff difference of ^ 40 K in the Balmer lines, is within the un- 
certainties of this analysis (including the errors from the obser- 
vations, model, broadening, NLTE effects, etc.). Barklem et al. 
(2002, hereafter PB02) discuss in detail the several uncertainties 
associated with the calculation of hydrogen lines, to which one 
should also add the uncertainties from the choice of hydrogen 
colHsions (PB07). 

The Pa/3 line, at 1281.8 nm, is not available in the Brault & 
Neckel (1987) disk-centre atlas nor in the Kurucz (2005) atlas, 
which only extend to 1250.8 nm and 1098 nm, respectively. For 
this line we use only the Kurucz et al. (1984) flux atlas, which 
is based on the same raw data as Kurucz (2005), but covers a 
somewhat larger wavelength range. However, the normalisation 
of the Kurucz et al. (1984) atlas in the region around Pa/3 is prob- 
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Fig. 7. Flux ratios between NLTE and LTE line profiles, for 
the hydrogen lines considered and for the different models. 
Wavelength difference AA measured from the line core. Results 
for the PHOENIX NLTE model atmosphere (not shown) are in- 
distinguishable from the corresponding LTE model in this figure. 
For Hff and Hy6 the ratio for the line core is not shown, to empha- 
sise the NLTE effects in the wings, formed in the photosphere. 



lematic. To compensate for this we re-normalised this part of the 
atlas, finding a suitable continuum level from a polynomial fit to 
nearby continuum high points. 



5.3. Synthetic profiles 

The computation of the H line profiles has been done allow- 
ing for departures from LTE. We use a 20-level hydrogen model 
atom (19 H I levels plus continuum) based on the atom of PB07, 
which was adapted from Carlsson & Rutten (1992). The colli- 
sional cross-sections by PB07 are employed. These include in- 
elastic collisions with electrons and hydrogen atoms (using the 
cross sections of Soon 1992, assuming E = ficm), mutual neu- 
tralisation and Penning ionisation. To speed up the calculations, 
especially in the 3D case, we have neglected the bound-bound 
transitions starting with a lower level n >- 6, thus including 
only 80 bound-bound transitions (out of a total of 171 in the 
original model atom). Tests with ID models show that the ef- 
fects of removing these lines are negligible. 

Line profiles were computed using our LTE code and the 
MPI-version of the MULTI3D code (Leenaarts & Carlsson 
2009). To save computational time, full 3D NLTE line forma- 
tion was performed only on eight snapshots of the simulation. 
Using MULTI3D to compute the LTE and NLTE line profiles 
we obtain the wavelength-dependent NLTE/LTE ratio, averaged 
for these eight snapshots. The final line profiles are then obtained 
by using our LTE code to compute the 3D LTE line profiles for 
all the 90 snapshots in the simulation and then multiplying them 
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by the average NLTE/LTE ratio for each wavelength. The very 
small variation between the NLTE/LTE ratios for the eight snap- 
shots indicates that this procedure is a very good approximation. 
For the ID models the computational requirements are unimpor- 
tant, and they are not time-dependent, but for consistency we use 
the same procedure of the NLTE/LTE ratio for each model. 

We present the hydrogen line profiles for flux (disk-averaged 
intensity) and disk-centre intensity. After convergence of the 
level populations, the flux profiles have been computed using 
a total of 32 inclined rays (8 yu-angles, 4 i^-angles); a vertical 
ray was used for the disk-centre intensity profiles. A rotational 
velocity of vvot = 1.8 kms"' was used for the disk-integration. 
For the ID models a microturbulence of 1 .0km s ' and a macro- 
turbulence of 2.5 kms ' (consistent with values derived from a 
sample of Fe i lines) were used, although these choices have lit- 
tle effect on the H line wings. The hydrogen opacity is calculated 
using the HLINOP routine (Barklem & Piskunov 2003; PB07). 
This ensures a proper treatment of self-broadening (following 
Barklem et al. 2000) and Stark broadening (following Stehle & 
Hutcheon 1999). 

5.4. Results 

We present the H Une profile results in Figs. 8 and 9 for flux and 
disk-centre intensity respectively. To quantify the differences be- 
tween the observations and synthetic profiles ax^ approach was 
carried out in the following way. First, regions close to the line- 
wing continuum were identified in the observations. These re- 
gions are indicated in grey in Fig. 8 and were used for both disk- 
centre and flux profiles. The differences between observations 
and synthetic profiles were calculated in these regions. For a tan- 
gible quantification of these differences, we have estimated the 
change in T^ff that each model would need to best match the ob- 
servations of each line. This was achieved using several MARCS 
models with 5500 < T^g < 6000 (afl with logg = 4.44 and 
[Fe/H]=0.0), whose hydrogen line profiles were calculated and 
used to derive an intensity ratio I{A, Tgff)/I(A, T^g - 5777 K). 
This intensity ratio was then multiplied by each synthetic line 
profile, to obtain the approximate line profile of each model 
for an arbitrary Jeff- This approach is only an approximation, 
as the variation of the hydrogen lines with Teff will vary from 
model to model. Nevertheless, it is good enough for this pur- 
pose. Using an optimisation procedure we calculated the Tf-g 
that, for each set of observations, minimises the reduced x^, de- 
fined as l/N • Yi Uobs - ^modei)^ /c^, where A^ is the number of 
wavelength points minus the degrees of freedom (one, in this 
case) and cr^ the measurement error, which we assume to be 
constant. These results are shown in Table 2. Also shown is the 
reduced x^ for the 7^.^ adjusted line profiles, summed over all 
the observations, and normalised by the value for the 3D model. 
This gives a measure of the goodness of the fits. 

To help connect the differences between models and observa- 
tions with the temperature structure of the models we calculated 
the temperature gradient, defined as: 



l.OOF 



VT = 



dlogior' 



(1) 



with the optical depth r evaluated at 500 nm. For the disk-centre 
intensity profiles of Ha, H/3, and Pay, the temperature gradient 
was calculated for three regions, corresponding to the forma- 
tion regions of three wavelength points in the wings of the line 
profiles. The three wavelength points Ai, A2, and A^ are defined 
as the velocity shifts from the line cores of respectively -609, 
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Fig. 10. Effect of multiple blends on the wings of H/3. Two syn- 
thetic Hy6 flux line profiles of a ID MARCS model are shown: 
using only hydrogen {dashed blue line), and including 220 other 
atomic lines {solid yellow line). Results for MARCS models with 
Teff + 100 K and -100 K are also shown {lower and upper thin 
black lines, respectively). 
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-316, and -196 kms '. For each wavelength the contribution 
function is calculated, and Vr taken as the average of a small 
depth range around the peak of the contribution function. For 
the 3D model this is done on a column by column basis, and the 
final value for VT taken as the mean of all the ID columns in all 
snapshots, weighted by the continuum intensity of each column. 
The resulting values of Vr, plotted against line depression of 
the disk-centre profiles, are shown in Fig. 1 1 . (The wavelength 
points used are also shown as vertical lines in Fig. 9.) 

5.5. Discussion 

5.5.1. NLTE effects 

The NLTE effects on the hydrogen lines are quantified in Fig. 7, 
where the NLTE/LTE flux ratios are shown. The main conse- 
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Fig. 8. Normalised flux profiles for the H lines HyS, Ha, Pay, and Pa/3. Compared with the solar observations of Kurucz (2005). The 
exception is Pa^, where the Kurucz et al. (1984) atlas is used because the Kurucz (2005) atlas does not cover these wavelengths. For 
this line the continuum was re-normalised (see text). Synthetic profiles were computed in NLTE. The regions used in the;if^ analysis 
are indicated in grey. 



quence of the NLTE effects in the H lines is a deeper core, 
which is of less interest here since it is little sensitive to the pho- 
tospheric stratification but rather determined by chromospheric 
conditions. However, as noted by PB07 and shown in Fig. 7, 
there are non-negligible eff'ects on the wings of the Balmer lines, 
making the wings weaker compared to the LTE case, at least 
with our particular choice of H collisions for the NLTE calcula- 
tions (see discussion in Barklem 2007)). For the Paschen lines 
the NLTE eff'ects are much smaller, and they cause only a deeper 
core. 



Compared with the results of PB07, we find a weaker NLTE 
effect in the wings of the Balmer lines. Using a similar MARCS 
model, the same recipe for the collisional rates, and a similar 
code for NLTE radiative transfer (MULTI), PB07 finds a max- 
imum NLTE excess flux on the wings of Ha of about 2.5%, 
whereas we find around 1%. The origin of this difference is 
our inclusion of line-blanketing for photo-ionisation transitions. 
This additional source of opacity, in particular in the UV, leads to 
a decrease in the photo-ionisation rates, bringing the level pop- 
ulations closer to LTE. Line-blanketing was not included in the 
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Fig. 9. Normalised disk-centre intensity profiles for the H lines Hf5, Ha, and Pay. Compared with the solar observations of Brault 
& Neckel (1987). Synthetic profiles were computed in NLTE. The three vertical lines correspond to the wavelengths Ai, Ai, and /I3 
used to make Fig. 1 1 . 

Table 2. Estimated effective temperature differences between the models and observations. 
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(ATeft) denotes the simple mean of the seven AJeff values. 

The;^'^ value is the sum of the individual reduced ^^, normalised by the 3D model results. 



MULTI version 2.2 used by PB07, but is in MULTI version 2.3 
and in MULTI3D. If we switch off line -blanketing in MULTI3D, 
we obtain essentially the same results as PB07. 



5.5.2. Temperature gradient 

The results for the temperature gradient VT in Fig. 11 show 
a correlation between the line strength and VT. Typically, the 
higher Vr, the stronger the normalised line profiles as one would 
naively expect. In most cases the Holweger & Miiller, (3D), 
and PHOENIX models seem to fall on the same linear rela- 



11 



Pereira et al.: How realistic are solar model atmospheres? 



tion, while the 3D and MARCS models, showing a similar rela- 
tion, have a lower VT for similar line strengths. When compared 
with the Holweger & Muller, the higher Vr of the PHOENIX 
and MARCS models is consistent with their predictions for the 
wings being much stronger than those of the Holweger & Muller 
model. However, the VT alone is not enough to explain the dif- 
ferences between MARCS and PHOENIX models, with the lat- 
ter having usually a larger VT but a smaller Tetj correction. 

Between the PHOENIX models, the differences in the AFeff 
corrections are consistent with what is seen in VT: the NLTE 
model has a higher VT, and consequently also larger AT^ff cor- 
rections (positive or negative). Between the Balmer lines, the 
temperature gradient of the PHOENIX models is lower for H/3. 
This is a consequence of the abrupt change in temperature gra- 
dient that these models show at logjQ tsqo ~ 0.3 (see Fig. 1). H/3, 
being formed deeper, is formed mostly on the flatter side of this 
knee, whereas Ha is mostly formed on the steeper side of the 
knee. This helps explain why these models predict Ha to be too 
strong, while predicting H/5 to be too weak. The MARCS model, 
that does not show this 'knee', predicts both Balmer lines to be 
stronger than the observations. 

Perhaps the most interesting departure from the VT rela- 
tion with line strength is the case of the 3D and (3D) models. 
These models have a very similar VT, but the predictions of 
the (3D) model are of much weaker lines than the 3D model. 
For all the tests in the hydrogen lines, the (3D) is closer to the 
Holweger & Muller model. This indicates that the spatial and 
temporal variations of the 3D model are important to describe 
the shapes and strengths of the hydrogen lines due to the very 
large non-linearity in the line formation of these high excitation 
lines. Ludwig et al. (2009) reached similar conclusions in their 
analysis of H lines using 3D CO^BOLD models and stressed 
that the "3D eff'ect" in terms of Teff depends sensitively on the 
particular adopted mixing length parameters in the ID model at- 
mospheres. 

5.5.3. Comparison with observations 

Overall, no single model seems to reproduce the observations 
perfectly. All of the models tested require different T^ff correc- 
tions for different lines. This is particularly true for the Balmer 
lines, whose Tf.fi corrections often have opposite signs. Most 
models predict a too strong Ho- line, while at the same time H/3 
is not strong enough compared to observations. Given the higher 
number of blending features in H/3 one wonders if this effect is 
not caused by a continuum depression because of all the blends, 
not considered in the hydrogen-only synthetic profiles. To make 
sure the single-line approximation is valid, we performed spec- 
tral synthesis in this region for the ID MARCS model, including 
220 additional atomic lines in the calculations. Data for these 
lines were extracted from the vald database (Piskunov et al. 
1995; Kupka et al. 2000), selected as the strongest lines in this 
region. The results, shown in Fig. 10, indicate that the blends 
have a negligible effect in lowering the local continuum of the 
wings of HyS, validating our single-line approximation. 

The MARCS model predicts line profiles that are too strong 
in all cases, requiring a negative Teff correction. The Holweger & 
Muller model suffers from the opposite effect: its predictions are 
considerably weaker than the observations, except for Ha, when 
they seem to be just about right. Both PHOENIX models seem to 
behave similary to the MARCS in Ha, but then have mostly pos- 
itive Tf.ff corrections for the Paschen lines, a likely consequence 
of the variations in their temperature gradients in the deepest at- 
mospheric layers, as discussed before. The 3D model predicts 



the wings of Ha to be stronger and the Paschen line wings to 
be weaker than the observations, but its prediction for HyS agrees 
very well with the observations. 

The Balmer line results of the MARCS model can be com- 
pared with the LTE results of PB02. Although PB02 used an 
earlier MARCS model (Asplund et al. 1997), slightly differ- 
ent regions for;(f^ with the 1984 Kumcz flux atlas, and possi- 
bly slightly different input physics, it is nevertheless relevant to 
compare our results with theirs. Calculating the LTE MARCS 
Balmer profiles for the same flux atlas used by PB02, we derive 
a Teff from each Balmer line. For H/3 our determined T^ff is in 
good agreement with PB02, but for Ha our Teff is 75 K lower. 
While difficult to pinpoint an exact cause, this difference is likely 
to come from differences in the input physics used to calculate 
the line profiles. 

Except for the (3D) model, the synthetic profiles of the the- 
oretical models agree better for flux than disk-centre intensity. 
The biggest difference in fitted Teft for intensity/flux is found for 
Ha in particular for the MARCS model. This is a hint of short- 
comings in the description of the solar temperature profile in the 
deeper layers of these models, as the disk-centre intensity pro- 
files are formed deeper than the flux profiles. The Holweger & 
Muller model shows a different trend: its predictions generally 
agree better for the intensity profiles, and the difference between 
the corrections for flux and intensity is smaller It is reassuring 
to find that the 3D model gives the smallest variation between 
flux and intensity corrections, again an indication of a realistic 
temperature profile. 

For all the line profiles considered, the 3D and the PHOENIX 
LTE and NLTE models give the smallest variation for (AT^eff), 
the mean of the Tetf corrections. While a measure of the internal 
consistency of the models, (AFeff) is a simple approximation that 
does not take into account the shape of the line profiles. With the 
highest x^ values of the models tested, the predictions of the 
PHOENIX models do not reproduce the shape of the observed 
line profiles very well, while the 3D model performs the best 
also in this regard. The Holweger & Muller is the model with 
the highest (Areff ). The difference (3D)-3D is significant for the 
hydrogen lines: a {AT^jf) difference of about 80 K, and a worse 
X^ for the (3D) model. Overall, the 3D model stands out from the 
other models. While not perfectly describing the observations, 
in particular Ha and Pa^S, it gives the smallest variation in Teff 
between flux and intensity profiles, one of the smallest variations 
of Teff among the different lines, and its predicted line shapes 
agree very well with the observations as evidenced by the lowest 
X^ among the models. 



6. Continuum intensity distribution 

6.1. Context 

Another relevant diagnostic of our simulations is the continuum 
intensity distribution and contrast, A/rms, of the granulation. A 
comparison of these with observations will reveal how well the 
simulations capture the differences in radiative transfer between 
the up- and the downflows. 



6.2. Observations 

Because solar observations are made with instruments with a fi- 
nite resolution and subject to other effects such as straylight, it 
is difficult to ascertain what the solar intensity distribution and 
granulation contrast really is. All instrumental effects need to be 
carefully considered in order to make meaningful comparisons 
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Fig. 12. Continuum intensity distributions for the solar disk-centre for three wavelengths, as a function of the normalised continuum 
intensity for the 3D model (thin blue solid line), the 3D MHD model (thin red solid line), compared with observations (thick yellow 
line). Also shown is the prediction from the CO^BOLD 3D solar model (dotted line) taken from Fig. 5 of WR09. The intensity 
distributions for our 3D models were averaged over all the snapshots. 



with 3D simulations, but accurately compensating for all optical 
and straylight effects is a difficult task. 

Several studies have found the observed A/in,s to be lower 
than that of the 3D simulations even after consideration of at- 
mospheric and instrumental seeing effects (e.g. Uitenbroek et al. 
2007; Kiselman 2008). For the Hinode Solar Optical Telescope 
(SOT, Tsuneta et al. 2008), Danilovic et al. (2008) have char- 
acterised in detail the instrumental effects for the Spectro- 
Polarimeter (SP) instrument, and Wedemeyer-Bohm & Rouppe 
van der Voort (2009, hereafter WR09) have done the same for its 
Broad Filter Imager (BFI). Both of these studies find that when 
the instrumental degradation is carefully modelled, the contin- 
uum intensity distribution and A/j^s agree very well with the pre- 
dictions from 3D models. WR09 in particular test several types 
of 3D models and find an overall good agreement with the space- 
based observations of Hinode/SOT. 

For our comparison we use the Hinode/BFI observations 
employed by WR09. These observations were taken in three 
BFI channels in three wide-band filtergrams with central wave- 
lengths of 450.45 nm, 555.05 nm, and 668.40 nm. The FWHM 
of the filter transmission profiles are respectively 0.22 nm, 
0.27 nm, and 0.31 nm. The observations were obtained in the pe- 
riod between November 2006 and February 2008. For our com- 
parison, we use the values from WR09 for the observations de- 
convolved with the instrumental profile (in an attempt to cancel 
out the image degradation), and compare them with the raw pre- 
dictions from the simulations (no image degradation applied). 

6.3. Results and discussion 

We compare our 3D models with the observations and the 
CO^BOLD 3D model results from WR09. CO^BOLD is inde- 
pendent of the stagger-code we employed, using different nu- 
merical methods and atomic physics. Results for the disk-centre 
intensity distributions at three wavelengths are shown in Fig. 12. 
The A/ims values for the same wavelengths are given in Table 3. 
Here A/in,s was calculated in the same way as WR09, follow- 
ing their equation (1), and averaged in time for all the snapshots 
considered. 

We find that our 3D model reproduces the observations well. 
Its A/„ns is sUghtly higher than for the CO^BOLD model, but 
otherwise results from the two models are very close, which 
is encouraging. The observed contrast is slightly higher than 



predicted from our 3D model, probably because (as noted by 
WR09) the observations seem to span a wider range of intensi- 
ties and in particular have a more pronounced 'tail' at high inten- 
sities. The double peaked structure of the intensity distribution 
represents the brightness from the inter-granular lanes (highest 
peak) and the granules (lower peak). Compared to the observa- 
tions, the 3D model has a pronounced minimum in the distri- 
bution between the two peaks, which can be attributed to the 
lack of magnetic features (such as bright points, ribbons, and 
other structures). The 3D MHD model, on the other hand, has 
a distribution that is closer to the observed, with no local mini- 
mum between the two peaks. This is probably because the bright 
points and other structures tend to blur the sharp transition from 
granule to inter-granular lane, and populate the distribution with 
intensities below the peak attained in the hottest granules. The 
reduced contrast of the MHD model is consistent with its dif- 
ferent stratification and shallower temperature gradient (warmer 
upper layers), when compared to the 3D model with no magnetic 
fields. Nevertheless, the 3D MHD model still does not show the 
high-intensity tail seen in the observations. 



Table 3. Disk-centre A/rms for the deconvolved observations, our 
3D model and the CO^BOLD 3D model. 



A/n,. (%) 




A (nm) 






450.45 


555.05 


668.40 


Observations'" 


26.7+1.3 


19.4 ±1.4 


16.6 ±0.7 


3D model 


25.4 ±0.8 


18.6 ±0.6 


14.3 ± 0.5 


3D MHD model 


24.1 ±0.8 


17.7 ±0.6 


13.7 ±0.5 


CO^BOLD" 


25.0 ±0.1 


18.1 ±0.1 


13.8 ±0.1 



^ From WR09. 



7. Fe abundances and line shapes 

7.1. Context 

The detailed shapes of photospheric absorption lines carry cru- 
cial information about the atmospheric conditions. The strengths 
of weak spectral lines mainly reflect the temperature structure 
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in the line-forming regions, at least in the framework of LTE, 
while stronger lines become increasingly sensitive to the veloc- 
ity field due to desaturation effects. All observed spectral lines 
show asymmetries due to the presence of (anti-)correlations be- 
tween temperature and velocities in the up- and downflows. The 
warmer upflows lead to high continuum intensities and blue- 
shifted, strong line profiles due to the steep temperature gradi- 
ents while lines from downflows are red-shifted, weak and have 
low continuum intensities; the line strengths are normally heav- 
ily biased towards the upflows, also because of their typically 
larger area coverage. The spatially averaged profiles thus be- 
come skewed, resulting in lines with typically blue-shifted cores 
and C-shaped bisectors for stars like the Sun (e.g. Dravins 1982; 
Asplund et al. 2000a, and Fig. 15). Naturally, ID hydrostatic 
model atmospheres are unable to explain such line asymmetries 
and furthermore require the introduction of two additional free 
parameters: micro-turbulence to obtain realistic broadening of 
partly saturated lines and macro-turbulence to get reasonable 
line widths even if the detailed shape of observed lines can never 
be fully recovered. With a self-consistent convective velocity 
field, 3D models like those employed here do not need to invoke 
micro- or macro-turbulence to obtain an excellent agreement 
with observed line shapes. Realistic 3D simulations achieve this 
from first principles, without recourse to adjustable parameters 
(Asplund et al. 2000a). 

Fabbian et al. (2010, 2012) have recently investigated the 
impact of magnetic fields on Fe spectral line formation in the 
quiet Sun, in particular in terms of the inferred solar photo- 
spheric Fe abundance. Using a series of 3D MHD simulations 
of varying magnetic field strengths computed with the Stagger- 
code, they investigated the 3D LTE line formation of 28 Fe i lines 
and found noticeable diff'erences: lines with small or negligible 
Zeeman-broadening are still affected by the different mean tem- 
perature structures in 3D models with magnetic fields. For an 
average vertical field strength of 10 mT, they found that the de- 
rived Fe i-based abundance is x 0.05 dex higher than without 
magnetic fields, the exact effect depending on the particular line 
in question. If this holds true, it would mean that the solar chem- 
ical composition presented by Asplund et al. (2009) would have 
to be revisited given that it was determined using the same non- 
magnetic 3D model as we are studying herein. 



7.2. Observations 

We make use of the solar FTS disk-centre intensity atlas of 
Brault & Neckel (1987) to measure the observed line shifts and 
line bisectors for a sample of Fei and Feii lines. The FTS at- 
las is on an absolute radial velocity scale as the relative mo- 
tion between Sun and Earth is known precisely and has been 
corrected for the solar gravitational redshift of 633 ms"' (for 
light intercepted on Earth, Lindegren et al. 1999). The line list 
stems from Asplund et al. (2009) and is augmented with addi- 
tional lines from Asplund et al. (2000a) The necessary laboratory 
wavelengths to place the measured line shifts and bisectors on a 
velocity scale comes from Nave et al. (1994) for Fei and from 
Johansson (1998, private communication, see also Nave 2012) 
for Fe II. 



7.3. Results and discussion 

We performed 3D LTE radiative transfer calculations of Fe i and 
Fe II lines using the 3D hydrodynamical solar model atmosphere 
also employed by Asplund et al. (2009) as well as a 3D MHD 




Fig. 13. Iron abundances derived from Fe i lines (circles) and 
Fell lines (triangles), as a function of excitation potential Ck'exc) 
for the 3D model (blue, filled symbols) and the 3D MHD model 
(red, open symbols). The lines show linear fits to the abundance 
relation for Fe i lines, for the 3D model (blue solid) and 3D MHD 
model (red dashed). 



simulation with an average vertical magnetic field of 10 mT from 
Thaler et al. (in preparation); we have also ensured that the corre- 
sponding T simulation of Thaler et al. produces Fe lines indis- 
tinguishable from those of the Asplund et al. (2009) model. The 
theoretical disk-centre intensity profiles have been spatially and 
temporally averaged; the time sequence corresponds to 45 min 
of solar time with snapshots every minute. No micro- or macro- 
turbulent broadening entered the 3D line formation calculations 
but the resulting averaged line profiles were convolved with 
a Gaussian corresponding to the finite resolving power of the 
Brault & Neckel (1987) solar atlas. Polarisation and Zeeman- 
splitting have not been considered in the radiative transfer calcu- 
lations. 

The solar Fe abundance has been derived from each of our Fe 
lines using both the 3D hydrodynamical model and the 3D MHD 
model. For simplicity, we have adopted the equivalent widths 
given in Scott et al. (2013, in preparation), and, when not avail- 
able there, in Asplund et al. (2000b). The magnetic fields can im- 
pact the line strengths in two ways: directly via Zeeman-splitting 
and indirectly via the different atmospheric stratifications, espe- 
cially the temperature structure. We are not in a position to ex- 
actly quantify the former effect but the latter will dominate for all 
of our lines due to their small Lande factors and relatively short 
wavelengths Fabbian et al. (2012). Had we considered Zeeman 
splitting the derived Fe abundances would be < 0.02 dex smaller 
for the few lines with non-negligible Lande factors in our line 
sample for the 10 mT simulation (Fabbian et al. 2012). 

Fig. 13 presents a comparison of the Fe abundances from 
the two 3D models. Because the MHD model is slightly warmer 
in the higher atmospheric layers (Fig. 1), one would expect the 
inferred Fei abundance to be higher than without considera- 
tion of magnetic fields in the convection simulation. Also, the 
Fe II abundances should show very small differences given their 
greater formation depths. We find both of these to be true: the 
MHD mean Fe i abundance is 0.06 dex higher than without mag- 
netic field, in line with the findings of Fabbian et al. (2012), 
while there is little difference for Fe ii. However, as is clear from 
Fig. 13 there is distinct trend with excitation potential for the 3D 
LTE Fe I abundances with the MHD model, which is not present 
with the model of Asplund et al. (2009). It should be noted that 
this trend cannot be removed by departures from LTE because 
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Fig. 14. Upper panel: Observed (crosses) line shifts for a sam- 
ple of Fe I and Fe ii lines for disk-centre intensity as a function of 
line strength together with predictions from the 3D model (blue, 
filled circles) and the 3D MHD ( 1 00 Gauss) model (red, open cir- 
cles). Lower panel: Differences between predicted and observed 
line shifts. 



those are positive and more so for low excitation lines, as demon- 
strated using the temporally and spatially averaged (3D) model 
(Bergemann et al. 2012; Lind et al. 2012). 

For both the observed and predicted line profiles, the line 
centres were determined using a cubic spline around the wave- 
length points with minimum intensity. The solar gravitational 
redshift (of light intercepted at Earth) of 633 ms"' was sub- 
tracted to obtain the observed central wavelengths on an abso- 
lute wavelength scale. Fig. 14 shows the observed and predicted 
Fe I and Fe ii line shifts relative to the adopted laboratory wave- 
lengths of the lines. In line with previous findings, weaker lines 
have a more pronounced convective blue-shift due to the larger 
depths of formation where convection and the anti-correlations 
between temperature and velocity are the largest; the cores of 
stronger lines become progressively less blue-shifted such that 
Fe lines with an equivalent width of ~ 10 pm have nearly van- 
ishing line shifts (Asplund et al. 2000a). The agreement between 
predicted and observed line shifts is very satisfactory for the 3D 
hydrodynamical model as demonstrated in the lower panel of 
Fig. 14: 30 + 60 ms"' for our Fei lines and -50 ± 70 ms"' for 
Fen. As also found by Asplund et al. (2000a) the stronger Fe 
lines tend to have slightly underestimated convective blue-shifts; 
Fe I lines with equivalent widths < 6 pm have a mean difference 
of only 9ms"'. Given the slightly deviating behaviour of two of 
the weakest lines (Fei 669.9 nm and Feii 562.7 nm) one could 
suspect that they are more affected by blends or erroneous lab- 
oratory wavelengths than the average line. The predictions from 
the 3D MHD model have the correct qualitative behaviour but 
have systematically too little convective blue-shifts; the mean 
difference for Fe i is 100 ± 50 m s"' . 

A comparison between observed and predicted line bisec- 
tors tell a similar story as the line shifts. Solar disk-centre inten- 
sity line profiles show a characteristic C-shaped bisector (weaker 
lines tend to show only the upper part) with a typical velocity 
span of 300-600 m s"' with the exact shape depending on the line 
formation height and temperature/velocity sensitivity (Asplund 
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Fig. 15. Upper panel: Observed disk-centre intensity bisectors 
for a sample of Fei and Feii lines. Lower panel: Differences 
between predicted and observed bisectors. The thick lines repre- 
sent the average difference over all bisectors. 

et al. 2000a). Fig. 15 shows the differences between the predicted 
and observed bisectors for our sample of Fe lines; ideally these 
differences should manifest themselves as vertical lines at zero 
velocity offset. The agreement is very satisfactory for the 3D hy- 
drodynamical model while the bisectors based on the 3D MHD 
are not sufficiently blue-shifted, in line with the Une centre com- 
parison. 

8. Conclusions 

Realistic solar atmospheres are of paramount importance for our 
understanding of not just the Sun but also of observations of 
other stars. The Sun provides an ideal test bench to test the phys- 
ical ingredients of the models, which if successful can then be 
applied to other stars with some confidence. A critical require- 
ment for a realistic model is that its thermodynamical quantities 
such as temperature, density and pressure match those of the real 
Sun. In this work we have undertaken a systematic study of the 
temperature structure of several solar models, using several key 
observational tests: continuum centre-to-limb variation, absolute 
continuum fluxes, wings of hydrogen lines, and also the inten- 
sity fluctuations over the granulation and detailed line shapes 
and asymmetries. 

In all diagnostics we find that the 3D model reproduces the 
observations very well. This is especially true for the centre-to- 
limb variations, where its remarkable agreement surpasses even 
that of the semi-empirical Holweger & Miiller model, which was 
built to fit the centre-to-limb variations. The 3D model also per- 
forms very favourably against the absolute continuum fluxes ob- 
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servations. For the hydrogen lines, the 3D model predicts the 
wings of the Ha line to be slightly stronger than the observa- 
tions, but on the other hand provides a very good agreement for 
the other lines, and the best overall agreement of all the models 
tested. In terms of the continuum intensity fluctuations over the 
solar granulation, it is reassuring to find that the 3D model repro- 
duces the observed intensity distribution and A/„ns well. The 3D 
model also predicts line shifts and asymmetries that agree very 
well with observations, which further supports its high degree of 
realism given the great sensitivity of the exact line shapes on the 
atmospheric conditions and line formation process. 

In light of the work of Fabbian et al. (2010, 2012), we also 
calculated the predictions of a simulation with an average verti- 
cal magnetic field of 10 mT (the 3D MHD model). Regarding the 
Fe line asymmetries, shifts, and abundances, the 3D MHD model 
agrees slightly less well with observations, suggesting that either 
the eff'ects of magnetic fields have been overestimated or that 
it is missing some ingredient that counteracts the consequences 
of the magnetic fields for the Fe line formation. Together with 
the evidence from the other diagnostics, it implies that at this 
stage there is no justification to prefer the solar abundances de- 
rived from the current generation of 3D MHD solar models over 
the 3D-based analysis of Asplund et al. (2009); our results sug- 
gest that the 3D MHD Fe abundance corrections advocated by 
Fabbian et al. (2010, 2012) are over-estimated. 

The ID theoretical models agree well with the observed 
absolute continuum fluxes, especially the MARCS model. 
However, both the MARCS and the PHOENIX models predic- 
tions for the centre-to-limb variations are consistently below the 
observations, both in the visible and in the infrared, which we at- 
tribute to a too steep temperature gradient. Such ID hydrostatic 
models obviously cannot predict any line asymmetries or inten- 
sity contrasts. We find that the small difference in the temper- 
ature structure between the PHOENIX LTE and NLTE models 
does not translate into any significant difference in our compari- 
son. Their results are very similar. If anything, the NLTE model 
performs slightly worse against the observational tests. This is 
likely to result from its somewhat steeper temperature gradient, 
due to NLTE cooling of the outer layers. 

The agreement between the predictions from the 3D model 
and the observations demonstrates its very high degree of real- 
ism in its temperature stratification. Together with its realistic 
velocity fields and treatment of convection as exemplified by the 
line asymmetries, it places the 3D modelling in an excellent po- 
sition to perform chemical abundance studies. It is noteworthy 
that it greatly outperforms any of the investigated ID models, 
both theoretical flavours such as MARCS and PHOENIX and 
the semi-empirical Holweger & Miiller model. There is thus no 
justification to continue to rely on the inferred solar abundances 
from ID-based analyses when a significantly improved alterna- 
tive is available. 
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